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Fracture is a fundamental mechanism of materials failure. Propagating cracks can exhibit a 
rich dynamical behavior controlled by a subtle interplay between microscopic failure processes 
in the crack tip region and macroscopic elasticity. We review recent approaches to understand 
crack dynamics using the phase field method. This method, developed originally for phase 
transformations, has the well-known advantage of avoiding explicit front tracking by making 
material interfaces spatially diffuse. In a fracture context, this method is able to capture both 
the short-scale physics of failure and macroscopic linear elasticity within a self-consistent 
set of equations that can be simulated on experimentally relevant length and time scales. 
We discuss the relevance of different models, which stem from continuum field descriptions 
of brittle materials and crystals, to address questions concerning crack path selection and 
branching instabilities, as well as models that are based on mesoscale concepts for crack tip 
scale selection. Open questions which may be addressed using phase field models of fracture 
are summarized. 



1. Introduction 

The dynamics of crack propagation is an important and long standing challenge in 
materials science and solid-state physics [TJ [2] , and in the recent years the physics 
community saw a rebirth of interest in the problem of dynamic fracture, also in 
combination with the concept of phase field modeling |3-7j. In this respect, fracture 
is yet another extension of a concept that was originally introduced to describe 
and simulate the kinetics of phase transitions - as the historical name "phase 
field" suggests. Today, this methodology has entered an "inflationary stage", and 
it has been successfully extended to various new applications in physics, materials 
science, but also biology and medicine. The characteristic property of phase field 
methods is the existence of one or more "phases" , for example solid or liquid phases, 
which can be transformed into each other. While a scalar phase field was originally 
introduced to distinguish between different phases O H], other fields have been 
introduced to characterize other properties such as the local crystal orientation [8]. 
On this broader scope, the phase field variables are therefore more appropriately 
denoted as order parameters. For the modeling of cracks, such an order parameter 
distinguishes between the solid phase and the "broken" state inside the crack. As 
usual in the phase field context, the order parameter changes smoothly between 
the states at the crack surfaces. The growth of a crack becomes then conceptually 
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comparable to the front propagation in a first order phase transition. During the 
past years, several important questions concerning the growth of cracks have been 
successfully tackled by this new paradigm. 

The uniform motion of a crack is relatively well understood in the framework of 
continuum theories [9l411|. Here, the conventional approach is to treat the crack 
as a front or interface separating broken and unbroken regions of the material; 
propagation is governed by the balance of the elastic forces in the materials and 
cohesive stresses near the crack tip |12H14| . Many characteristic features of crack 
propagation are nowadays well established by experimental studies [T5H22] . As soon 
as the flux of energy to the crack tip exceeds a critical value, the crack becomes 
unstable and starts to branch while emitting sound waves. These phenomena are 
consistent with the continuum theory of sharp crack tips, but it fails to describe it, 
because the details of crack growth, in particular the chosen crack path and veloc- 
ity, depend on details of cohesion on microscopic scales |23j . Nevertheless, empirical 
energy balances and simple propagation laws that are frequently used in engineer- 
ing applications, cannot account for the richness of actual fracture phenomena. In 
particular, they cannot predict dynamical instabilities of fast moving cracks. The 
fundamental mechanisms of those instabilities have been extremely difficult to elu- 
cidate because they appear to result from a non-trivial coupling between dynamical 
phenomena inside the crack tip region, known as process zone, and linear elasticity, 
with no clear separation of scale between atoms and the system size. 

Large scale molecular dynamics (MD) simulations with about 10 7 atoms allowed 
to get deeper insights into the growth behavior of cracks [24-281. Although limited 
to submicron samples and very short timescales, these simulations were able to 
reproduce key features of crack propagation like the initial acceleration and the 
onset of instabilities. Nevertheless, a detailed understanding of the complex physics 
of crack propagation, in particular aspects of the pattern formation process, still 
remain a major challenge |29j . 

At this level, continuum descriptions, in particular phase field methods that 
avoid dynamical artifacts which are associated with the breaking of translational 
and rotational symmetry fESEX]; offer a useful and complementary perspective on 
crack propagation as a pattern formation process. The focus of this article is on 
phase field formulations. Fracture has also been modeled by level set approaches, 
see e.g. [32]. Furthermore, cracks are only a particular type of defects in materials; 
for a recent overview on the phase field modeling of defects in general see (33] . 



2. Basic aspects of fracture 

One of the cornerstones of fracture in brittle materials is the Griffith criterion 
|34| . Materials under tension store elastic energy with density w ~ P 2 /E with 
the applied stress P and an elastic modulus E. The appearance of a crack with 
length L leads to a local relaxation of the elastic strain, AF e i ~ wL 2 (here in two 
dimensions), but at the same time the new crack surfaces are created, which lead 
to an additional surface, or more generally, fracture energy, AF S ~ 7L, see Fig. [TJ 
For short cracks, the increase of the interfacial term is most important, whereas 
for long cracks the second term dominates; hence fracture is similar to nucleation 
in first order phase transitions, eventually with a cohesive interaction between the 
crack surfaces [35]. The critical length of a crack, Lq ~ E^f/P 2 , is called the Griffith 
length. 

Near the crack tip the stress has a universal singular behavior in the framework 
of the linear theory of elasticity [361 EZ] > see Fig. [lj 
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Figure 1. Stress relaxation and interface creation due to the appearance of a crack. Near the (sharp) 
crack tip the stress becomes singular in the vicinity of the crack in the framework of the linear theory of 
elasticity. 




where K m are the stress intensity factors for the three modes of loading (see Fig. [2]), 
and 9 is the angle between the radial vector of length r with origin at the (sharp) 
crack tip and the local crack direction. The angular depedences f™{&) are known 

functions, see e.g. [9rlll|. We note that the strength of this singularity, a ~ r -1 / 2 , 
differs from dislocations, where stresses scale as a ~ r . In both cases, the cutoff 
for the singularity and the regularization of stresses is important since stresses 
cannot be infinite |38j ; this can be due to nonlinear phenomena in the process zone 
or a finite crack tip radius. 

For propagation of a crack, the energy release rate, also known as crack extension 
force, 

G= 1 -^(K] + Kh) + ^K] II (2) 

has to exceed a material dependent threshold G c ; this is equivalent to the Griffith 
criterion. In the ideal case, this energy barrier is exactly twice the surface energy 
of the material, G c = 27, since two new surfaces are created; in reality however, 
this threshold is usually much higher. 
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3. Continuum field models 

The main aim of continuum field models of fracture is to provide a unified frame- 
work for the whole phenomenology of fracture, ranging from crack initiation to 
oscillations and branching. 

Probably the first model in a series of phase field like descriptions of crack growth 
in brittle materials was developed by Aranson, Kalatsky and Vinokur |39j . focus- 
ing on mode I cracks in two dimensions. Apart from the displacement field, which 
obeys the elastodynamic equations together with a damping term, they introduce 
an order parameter p which accounts for the dynamics of defects. This order pa- 
rameter has the value p = 1 in the solid, whereas it vanishes inside the cracks, 
p = 0. It is assumed that this order parameter changes smoothly at the crack 
surfaces on a scale that is large in comparison to the atomic spacing, to justify a 
continuum description. This phenomenological model already captures many im- 
portant features, like crack initiation, propagation, dynamic fracture instability, 
sound emission, crack branching and fragmentation. 

In this model, the elastic modulus is assumed to be proportional to the order 
parameter, E = Eqp, which implies that the interior of an ideal crack is stress free. 
The stress is given by 



where the last term mimics a hydrostatic pressure due to the creation of new 
defects. The elastodynamic equation contains a viscous damping 



The evolution of the order parameter is assumed to follow a relaxation law 



which contains several material parameters. The first, diffusive term is analogous 
to the contribution of the gradient term in a phase field model, and similarly the 
second term corresponds to a double well potential with minima at p = and p = 1. 
The last term reflects an advective contribution, which couples the local order 
parameter to the displacement rate and is responsible for the localized shrinkage 
of the crack due to material motion. It turns out that this phenomenological term 
is important to maintain a sharp crack tip. 

For low driving forces, the model predicts steady state growth of a single straight 
crack in a strip, and the distribution of the elastic fields is in agreement with linear 
elastic fracture mechanics. Above a critical crack velocity, which is some fraction 
of the Rayleigh speed vr (sound speed of a wave traveling along a free surface), 
branching occurs; the threshold velocity depends on the material parameters and is 
roughly in the range (0.3 — 0.55)t>ij. The instability manifests itself as pronounced 
velocity oscillations, crack branching and sound emission from the crack tip. A 
velocity gap, i.e. a minimum crack velocity as in lattice models |23| . is not found. 

A frequent and important analysis for phase field-like models for crack growth 
in a channel geometry is related to the asymptotic behavior in the tail region 
far behind the crack tip. There (depending of course on the boundary conditions, 
and we assume that the strip is loaded by a fixed displacement), both the order 





p Ui = 7]V 2 Ui + 



d(Jjj 



(4) 
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Figure 3. Steady and unsteady propagation of mode III cracks with increasing load from (a) to (c). 
Straight crack propagation is stable at low velocity (a), becomes oscillatory at intermediate velocity (b), 
and exhibits branching (c) at high velocity. Taken from I41| . 



parameter and the elastic fields become homogeneous, and therefore a theoreti- 
cal analysis becomes feasible. In the above model |39| the crack opening depends 
logarithmically on the strip width L for fixed applied strain. This is in contrast 
to the expectation that the material should be fully relaxed in the wake, thus all 
displacement should be accumulated within in crack, and hence the crack opening 
should increase linearly with the strip width. 

To overcome this problem, Karma, Kessler and Levine (KKL) proposed a differ- 
ent phase field approach where the scalar order parameter describes the state of the 
material in Lagrangian "material coordinates", instead of the density in Eulerian 
coordinates [ID]. They investigated in particular mode III cracks, which makes the 
description simpler since only the out-of-plane displacement component u z exists 
in a two-dimensional description. They link the description of cracks closer to con- 
ventional phase field modeling, as it is widely used e.g. for solidification. The order 
parameter distinguishes between broken, = 0, and unbroken states, <p = 1, and 
its dynamical evolution is derived variationally from a free energy functional, 



dx 



(6) 



and rdt4> = —5F/5(p. The "bulk" states correspond to the minima of a double 
well potential, Vdw = </> 2 (l — 0) 2 /4, which together with the gradient term is 
responsible for the fracture energy. The elastic term involves the shear modulus 
[i and the strain e, which induces breaking of the material if it exceeds a critical 
strain e c . The choice of the coupling function g((j)) turns out to be critical and 
should scale as g((j)) ~ cj) 2+a with a > for small values of cj), in order to obtain 
full stress relaxation in a completely broken solid in the limit of large system 
sizes. We mention the similarity of the model to conventional phase field models of 
solidification, where the last term corresponds to the deviation from the melting 
temperature g(4>)L(T — Tm)/Tm', this also shows that the crack motion is here a 
fully reversible process, i.e. the crack can heal if the load is released. The choice of 
the coupling function, however, is here more strict, which is due to the fact that 
the crack width is directly related to the phase field interface thickness, hence the 
model does not operate in or close to the sharp interface limit. 

The model was applied to the dynamics of mode III cracks in [5T] . If the relax- 
ation timescale r is long, the growth of the crack is slow and controlled by the 
interface kinetics; in this limit the elastic fields are quasistatic. In contrast, for 
fast growth the crack speed is limited by the sound speed, and all dissipation in 
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this brittle limit takes place in the interface region. Three different basic regimes 
of solutions are found, see Fig. |3| Steady state growth for low driving forces, an 
asymmetric tip-splitting mode for intermediate loads, where the crack tip follows 
a snake-like sinusoidal trajectory, and finally a chaotic tip-splitting regime with 
well-developed branches for large loads. The maximum velocity of a single straight 
crack in the inertial limit turns out to be v max « 0.41c s (c s is the shear wave speed), 
and it is argued that the limitation of the crack speed is due to tip blunting as a 
result of relativistic contraction; the observed branching angles close to the onset 
of branching are similar to analytical predictions [42J . 

Henry and Levine extended the model further to two-dimensional plane strain 
situations and investigated the growth of cracks under mode I and mode II condi- 
tions [43] . recovering the same basic dynamical behavior of the previous mode III 
phase field study of fast moving cracks. Also, they generalize the elastic energy and 
also suggest a symmetry breaking term to go beyond situations where the elastic 
energy is quadratic in strain, where the situation is symmetric under tension and 
compression without advective contributions, implying that cracks also grow in 
compressed materials. Therefore, Henry and Levine generalize the contribution of 
the elastic energy in the free energy functional, that is used to express the dissi- 
pative phase field dynamics, g(<p)(E ( f > — Eq) (Eg corresponds to in the mode III 
model), to 

w _/ I A 4 + A*4- iftre>0, 

with the modulus of compression K = (A + /x)/2 and a > 1 is a parameter that pre- 
vents breaking under compression. This avoids situations where additional shear 
forces on top of a mode I loading lead to splitting of the crack and symmetrical 
growth of the two branches around the previous propagation direction. Instead, the 
material breaks now only in the tensile region. For higher loading also the mode I 
crack becomes unstable, and secondary cracks develop out of a branching instabil- 
ity. Notice, that not only the symmetry breaking model by Henry and Levine but 
also the original KKL model satisfy the principle of local symmetry, i.e. a reorien- 
tation of the crack growth direction such that the mode II stress intensity factor 
vanishes, Kjj = 0. This property will be discussed in more detail in section [4} 

The model was then applied to biaxial loading (an additional tensile stress com- 
ponent acts in the growth direction), where a supercritical Hopf bifurcation from 
straight to sinusoidal oscillatory is observed. A similar phenomenon was previously 
observed experimentally in rubber [44]. Remarkably, this instability, which may be 
interpreted as the attempt of the system to get rid of the additional elastic energy 
related to the longitudinal strain, is observed only for dynamical crack growth but 
not for overdamped elasticity; the wavelength of the instability depends linearly 
on the strip width that is used for the simulation. 

A recent and more detailed investigation of the model, especially concerning the 
onset of branching, has been done in [35] , and a reasonable agreement of critical 
velocities and branching angles in experiments [THl 13 0B] and theoretical predic- 
tions |47[ I4"B] has been found. Also, the dynamics are altered in comparison to the 
preceding models such that dt4> < 0, which means that cracks can only grow and 
not recede. Therefore, sidebranches do not disappear even far behind the crack tip, 
otherwise they would shrink in order to minimize the interfacial energy. The same 
concept, together with the use of nonlinear elastic laws, is briefly presented in [49J. 

Marconi and Jagla presented a diffuse interface approach to slow fracture in brit- 
tle materials using a different concept [50]. They do not introduce an additional 
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order parameter to discriminate between broken and unbroken regions, but use the 
strain field itself as "order parameter" . The description of material failure is here 
entirely through the nonlinear form of the elastic free energy density, which satu- 
rates for high strains, therefore representing the broken material inside the crack; 
the limiting value of the nonlinear stress-strain relation determines the fracture 
energy. A similarity to phase field models stems from the presence of a gradient 
energy term, which acts here on the components of the strain tensor, thus intro- 
ducing also second order derivatives of the strain in the equations of motion (in 
other models, the strain appears in the phase field evolution equation only locally 
and with first derivative in the elastic equation). Within this model, the fracture 
energy depends on the load, and ad hoc modifications of the model are necessary 
for corrections, which introduce additional parameters to the theory. 

Most phase field models for fracture use nonconserved order parameters, and 
exceptions are rare. One of them is the work by Eastgate et al. [51]. In their case, 
the phase field is interpreted as the (normalized) mass density of the material 
surrounding the crack, similar to |39j . Inside the crack, it becomes zero, whereas 
in the solid it is = 1 — e^, which takes into account the local compression 
or expansion of the material due to elastic strain. Both the dynamics of the phase 
field and the elastic displacement are derived from a free energy functional, and the 
order parameter evolves according to 4> = — V-j, where the flux j is not only driven 
by the energy decay, but also contains the advective contribution due to the elastic 
displacement. The displacement field obeys a viscous law, and is therefore mostly 
appropriate for the description of fracture e.g. in colloidal crystals. This is one of 
the main differences in comparison to the work of Bhate et al. [52] for modeling 
of stress voiding in electromigration, where quasistatic elasticity is assumed. The 
model |51j is used to describe the deformation of an initially spherical hole in a 
two-dimensional strained system and the evolution towards a crack. Limitations 
of the model are related to the fact that inside the crack the phase field does not 
fully reach the vacuum state = 0, but retains a value that is proportional to 
the tangential strain along the crack surfaces. Related to this artifact is a slight 
deviation from the correct value of the Griffith point. 

4. Crack path prediction 

The prediction of the path that a crack chooses while it propagates through a 
brittle material has been a long standing problem in fracture mechanics. In a 
conventional picture, the equations of linear elasticity are solved for traction free 
crack surfaces with sharp crack tips [10] • Although the Griffith condition provides a 
criterion for crack growth, it is insufficient to predict the curvilinear crack paths or 
crack kinking, or even branching angles. The generally accepted condition Ku = 
assumes that the crack propagates in such a direction that it is in a pure opening- 
mode state with a symmetrical stress distribution about its local axis [53]. This 
principle of local symmetry has been rationalized using plausible arguments in 
the classical framework of fracture mechanics |54l [55] and has also been argued 
from the continuity of the chemical potentials at the fracture tip [56 1. Variational 
formulations based on energy minimization were considered in |57|, I5B] as a proposal 
to overcome the limitations of Griffith's theory. Hodgdon and Sethna suggested 
equations of motion for a crack based on symmetry arguments |59| . Nevertheless, 
crack propagation laws are hard to derive without a detailed knowledge of the 
physics of the process zone, where the elastic energy is transformed into new crack 
surfaces and dissipated. In addition, symmetry considerations cannot be invoked 
when the fracture energy and/or elastic properties become anisotropic |60j . 
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Hakim and Karma used the KKL phase field model |40j to derive the principle of 
local symmetry in a fairly generic way, and provide even generalized conditions to 
predict crack paths [6"2~] . Similar to considerations in classical fracture mechanics 
they calculate the energy flux into the crack tip. Apart from a dissipative term, 
it consists of purely elastic contributions, which are conservative and can therefore 
be expressed by the energy flux through a circle around the crack tip; since its 
radius can be large, this contribution to the energy - excluding the contribution 
from the segment that crosses the crack surfaces - is described by the singular stress 
field and yields the same result as Rice's J integral [63J, which is the crack extension 
force G, and also Eshelby's torque Gg = dG{6)/d9, where G{9) is the extension 
force for a slightly kinked crack into direction 9 |64j . This torque term tends to 
turn the crack in a direction that maximizes the crack extension force. The crucial 
difference to classical fracture mechanics is that the phase field model also contains 
a description of the process zone, and the integral contribution along a line that 
crosses the crack gives contributions from cohesive forces, thus being related to the 
crack surface energy. A similar torque term gives a contribution proportional to the 
Herring torque 70 = d'y/dO [65], and therefore turns the crack into a direction that 
minimizes the surface energy in anisotropic media. The main result for quasistatic 
fracture is the generalized principle of local equilibrium, stating 



where f2 is the dissipative force perpendicular to the crack, which vanishes iden- 
tically in isotropic media and otherwise becomes negligible close to the Griffith 
point. For isotropic media, 70 = 0, the prediction recovers therefore the antici- 
pated relation Ku = 0, otherwise a shear loading balances the Herring torque 
term. 

The principle of local symmetry and an alternative suggestion, the principle of 
maximum energy release rate, give identical results for smooth curvilinear cracks 
and only small differences for kinks [541 [66l [67] , and thus it is difficult the discrim- 
inate between the predictions. In anisotropic media, however, these principles can 
give substantially different predictions for the kink angles for specific surface energy 
anisotropy and loading conditions. The phase field simulations in [62] show that 
the kink cracks emerge from the main crack tip at an angle, which is initially close 
to the prediction of the maximum energy release rate criterion, but approaches 
the angle predicted by the force balance condition (generalized principle of local 
symmetry) on the scale of the process zone, see Fig. |4j 

For antiplane shear, where only one stress intensity factor exists, the prediction 
of the crack path requires consideration not only of the singular stress term, but 
also of the following, subdominant contribution in the stress field expansion, and 
the principle of local symmetry demands here that its prefactor becomes zero [68J . 

The problem of crack path prediction is also essential for the case of thermal 
loading, and has been investigated intensively theoretically and experimentally. 
Typically, a glass strip with a notch at its end is pulled at a constant velocity from 
an oven into a cold bath. The control parameters in this experiment are mainly 
the width of the strip, the temperature gradient between the oven and the cold 
bath, and the pulling velocity. If the velocity is small enough, a crack does not 
propagate. Above a first critical velocity, the crack starts propagating following a 
straight centered path, and above a second critical velocity, the crack begins to 
oscillate with a well defined wavelength. For a recent investigation, including a 
phase field model for this problem, which shows that the transition from straight 
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Figure 4. Comparison of phase field simulations and theoretical predictions of the maximum energy release 
rate criterion and the force balance condition derived from the phase field model for kinked cracks in a 
material with an anisotropic surface energy. Taken from |62| . 



to oscillatory cracks is supercritical, we refer to [69] and references therein. 



5. Phase field modeling of fracture in crystalline materials 

Crack propagation in crystalline materials is inherently anisotropic because the 
fracture energy depends on the orientation of the fracture surface with respect to 
some underlying set of crystal axes. For this reason, cracks are often observed to 
propagate along low energy cleavage planes. Crystalline anisotropy can in principle 
be included in the KKL phase field model [40J by modifying square gradient terms 
in the free-energy density to break rotational invariance |6H 162] . Using the general 
propagation law for crack paths in anisotropic materials given by Eq. ([8]), Hakim 
and Karma have derived a condition for a crack to escape from a cleavage plane. 
This condition stems from the physical picture that a crack will be "trapped" along 
a low energy plane corresponding to a cusp in the 7-plot (i.e. the plot of surface 
energy as a function of orientation) until the Eshelby configurational torque and 
hence Kjj is large enough to rotate the crack away from this plane. 

Another independent development of crack propagation in a phase field context 
has been developed by Jin, Wang and Khachaturyan [70] in the spirit of the phase 
field microelasticity approach, that had previously been used to study e.g. coher- 
ent precipitation and ordering [71H73] . martensitic transformations \74\ 175] and 
dislocation dynamics |76l 177] . The main step is here the solution of the elasto- 
static problem, since the discontinuous system, that contains cracks and voids, is 
not elastically homogeneous. Instead the elastic modulus of the material becomes 
equal to zero within the inclusions. For the solution of this problem beyond a per- 
turbative treatment [78J, the material is equivalently considered as a continuous 
homogeneous body with a mesoscopically heterogeneous misfit-generating stress- 
free strain e?-, also know as eigenstrain. If the misfit strain is known everywhere, 
the elastic problem is solved by minimizing the free energy E e i with respect to 
the displacements. Here, the eigenstrain is an additional degree of freedom that is 
present only inside the cracks or voids. Additional minimization with respect to 
the eigenstrain leads to a stress- free state inside the inclusions [79 1. This is con- 
ceptually similar to the representation of a crack as pile-up of dislocations [80J, 
and the eigenstrain e^-(r) plays the role of a long-lange order parameter. It can be 
determined by a time-dependent Ginzburg-Landau equation, which is first applied 
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inside the cracks and voids only, 

dt ~ LijU 5e%{v,t) (9j 

with positive definite Onsager coefficients Lyjfcj. 

To include also a dynamical evolution of the cracks, it is assumed that the cracks 
can develop along cleavage planes a = 1 . . .p with normal vectors H(a). The stress- 
free transformation strain is then assumed to have the structure 



e° lj (r) = ^h l (a,r)H J (a), (10) 



where the hi play now the role of the phase fields. An additional local energy term, 
which depends on the order parameters hi represents a cohesive energy between 
the crack lips and can be tuned to the desired model. Similar to conventional phase 
field models, a gradient energy term is added, 



P r 

E gra d = y2 /<?rod(h(a,r))(fr, (11) 
a=l Jv 



and the gradient energy density reads in the simplest case that the crack tip energy 
does not depend on the direction of the crack front in the cleavage plane, 

fgrad = /%)[H(a) x V(h(a,r) • H(a))] 2 (12) 

with material parameters /3(a). The evolution of the long-range order parameters 
follows then a usual relaxation equation from the total free energy F 

dhi(a,r,t) 5F 
at ohj(a,T,t) 

with an additional Langevin noise term £j. This shows the conceptual similarity 
to models for dislocation dyamics \76\ [77], where the role of the functions h is 
played by the Burgers vector. In contrast to the preceding phase field models, this 
description leads to crack paths that follow the cleavage planes in a crystal, and 
therefore do not obey the principle of local symmetry as in isotropic materials. 

The interactions of dislocations and martensites with free surfaces, voids and 
cracks have been demonstrated to work robustly in the framework of these models 
e.g. in |81H83j . An extension of this model to phase field simulations of crack 
tip domain switching in ferrorelectrics is presented in [84, [85] . These materials are 
widely used to fabricate sensors and transducers, but tend to break for high electric 
or mechanical loading due to their brittleness. Hydrate precipication and delayed 
hydrate cracking in zirconium are considered in [86J; here hydrogen or deuterium 
diffuse along a stress gradient towards a crack tip, where hydrides form and grow, 
and finally the crack develops through the hydride. 

Crack propagation has also been simulated using the phase field crystal model 
inspired from classical density functional theory, which resolves spatially the crys- 
talline density field [87J. Since the phase field crystal model naturally describes 
dislocations, it could potentially be used to investigate the transition from ductile 
to brittle crack propagation. However, further investigations will be necessary to 
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Figure 5. Crack propagation using the phase field crystal model. The greyscale illustrates the energy 
density and visualizes the stress relaxation in the vicinity of the crack. Taken from |87] . 

elucidate the applicability of the model for fracture, which remains almost com- 
pletely unexplored beyond pictorial examples of crack branching shown in Fig. [5] 
One limitation of the original phase field crystal model is the restriction to only 
one (diffusive) timescale, which limits not only the interface velocity but also the 
elastic or plastic relaxation in the bulk. This is appropriate for slowly propagat- 
ing cracks, but is restricted in applicability for higher velocities, since then elastic 
waves should travel with the speed of sound instead of relaxing diffusively. A po- 
tential method to decouple these processes is to use an acceleration term in the 
phase field crystal evolution equation, as demonstrated in |88j . 

6. Models with sharp interface limit 

Other phase field models of fracture have been developed from a different viewpoint 
than the models described so far in this article. Closer to the original motivation for 
introducing the phase field model in the context of phase transformations, this class 
of models uses the phase field method as a numerical tool to solve free-boundary 
problems for fracture. Those free-boundary problems have been formulated in an 
idealized picture of fracture, which is viewed as the highly nonlinear stage of the 
Asaro-Tiller-Grinfeld (ATG) instability of a stressed solid surface |89H91j . In this 
picture, the crack dynamics is governed by a well-defined set of sharp-interface 
equations. In contrast to the other models, the crack tip scale is not set by the 
phase field interface thickness, but by a "physical" selection principle that can be 
rigorously derived from the analysis of the sharp-interface equations. This approach 
has close analogies with the dendritic growth problem, where the selected tip radius 
is of the order of the geometric means of the capillary and the diffusion length. 
However, due to the different physics of the crack problem, the selection of the 
tip scale and growth velocity is different in details. For both dendrite and crack 
growth, the product of the tip radius and the growth velocity is fixed by the 
driving force (undercooling or applied elastic load, respectively), which determines 
the rate of energy dissipation. For dendritic growth, the steady state velocity scales 
proportionally to the ratio of heat or solute diffusivity and a microscopic capillary 
length, and is nontrivially selected through the anisotropy of surface tension. In 
contrast, for fracture, this velocity scales proportionally to the sound speed and 
is selected even for an isotropic fracture energy. Crack tip selection requires no 
adjustable parameters and naturally predicts propagation velocities appreciably 
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Figure 6. The Asaro-Tiller-Grinfeld instability. The (two-dimensional) isotropic solid is stretched or com- 
pressed with a uniaxial stress P along the interface, which leads to a morphological long-wave instability 
of the front. Material redistribution can be due to surface diffusion (left) or, if the solid is in contact with 
its melt, solidification and melting (right). 



below the Rayleigh speed. However, those predictions are generally limited by the 
strong assumptions made in formulating the sharp-interface equations governing 
the crack dynamics, and it is not yet clear how to relate those predictions to 
experiments for specific materials. 

The starting point for this approach is based on the ATG instability, which 
predicts that a uniaxially strained solid surface is morphologically unstable [89- 
|9T] , see Fig. [6j The reason is that a deformation of the surface by a rearrangement 
of material can lead to an overall reduction of the total free energy. Although 
the corrugation of the surface increases the surface free energy, the elastic energy 
can be reduced even more, since the stresses are released close to the interface. 
The dynamics of this process is usually assumed to be either surface diffusion 
for a free surface or melting-solidification dynamics for a solid in contact with 
its melt at or close to the melting temperature. Here we mention that an elastic 
deformation of the solid increases its free energy relative to the liquid, which can 
induce a stress-induced melting process. For a perturbation of the interface contour 
y(x, t) = A cos kx exp(Ai) the spectrum of the instability is for melting-solidification 
dynamics in linear approximation 

\ = K[kL ATG -(kL ATG )\ (14) 

with a kinetic coefficient K for the interface kinetics and the Grinfeld length 

Ej 

LATG= 2(1-^' ^ 

with the surface energy density 7, Young's modulus E and the Poisson ratio v. 
Notice that this lengthscale is - apart from dimensionless factors - the same as 
the Griffith length for crack growth, since both processes are the result of a com- 
petion between a release of elastic and an increase of surface (or fracture) energy. 
As a result, long- wave perturbations are unstable, whereas short-wave corrugations 
are suppressed due to capillary effects, which are represented by the second term 



in Eq. (14). The initial and late stage of the ATG instability has been modeled 
analytically, with sharp interface and phase field descriptions, see e.g. [92H97j and 
references therein; recently, also phase field models for surface diffusion reproduced 
the instability [981 [99] , and also phase field crystal investigations have been pre- 
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formed [TOO]. 

In the late stage of the instability, the corrugations form deep grooves, which 
accelerate and evolve toward crack-like notches in the solid, see Fig. [7j In the 
framework of the static theory of elasticity, the tips of these "cracks" become 
arbitrarily sharp and reach an infinite velocity if no other cutoff mechanisms are 
provided (finite time cusp singularity). This fast growth, apart from the missing 
regularization for static elasticity, therefore proposes a deeper link between the late 
stage of the ATG instability and crack growth. 

The first rigorous connection between cracks and the ATG instability was high- 
lighted in [56], where it was observed that the conditions for the occurrence of 
the instability, a uniaxial stress, are satisfied along the crack surfaces. For a mode 
I crack with finite length, that it subjected to a tensile stress far away, the nor- 
mal and shear stresses vanish on the crack surfaces, but the tangential stress is 
different from zero. Therefore, the appearance of the instability is predicted for 
sufficiently long cracks, since the condition of fixed crack tips leads to "quantiza- 
tion" conditions for the spectrum of the instability, and at least one unstable mode 
has to fit on the straight crack. The exact length threshold for the instability was 
later calculated numerically in [101J, which shows that this long- wave instability 
can be expected for cracks that are already several times longer than the Griffith 
threshold, i.e. for growing cracks. 

Up to this point, the dynamics of the crack tip has been excluded from the 
considerations, and in |102] a first model was proposed to describe the motion of the 
crack tip itself as a stress induced rearrangement of material. There, in particular, 
the advancement of the crack was suggested to follow effectively a surface diffusion 
process. The crack tip was introduced as a new degree of freedom in the sense that 
it has a spatially extended structure (finite tip radius), which removes the problem 
of the stress singularity. Crack propagation therefore requires the redistribution of 
material if such a crack penetrates a solid. If the excess material (the volume inside 
the crack) would have to be transported out of the crack completely, the process 
would be slow. For a short diffusion path of the order of the tip scale, however, 
this is a conceivable mechanism for fast and steady state crack propagation. This 
suggested growth mechanism, which is nothing else than the late stage of the ATG 
instability, therefore requires a scale selection of the tip, in order to prevent the 
aforementioned finite-time cusp singularity. It was suggested that for fast crack 
growth the limitation to the Rayleigh speed provides the selection of the velocity, 
and together with the relationship between driving force and dissipation rate (the 
analogy of the Ivantsov relation), also the tip scale is selected. Hence the only 
difference in comparison to the conventional treatment of the ATG instability is to 
take into account the finiteness of the sound speed, and this leads to a regularization 
of the theory. In [102J , the model was explored within the framework of a local stress 
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Figure 8. Steady state growth of a crack with finite tip radius, determined through the interface kinetic 
coefficient and the sound speed. Taken from [103] . 

model, which does cover potential tip instabilities. 

A complete modeling of the problem using a phase field model was done in |103j . 
In contrast to |102j . a phase-transition process is assumed at the crack surface, 
since it can be modeled more easily than the than the higher-order surface diffusion 
process [HH1 US] • The general selection principles, however, are analogous for both 
growth mechanisms. There it is concluded that the suggested selection mechanism 
(finite sound speed) should indeed allow steady state growth and therefore prevent 
the finite time cusp singularity, thus steady state solutions with finite crack tip 
radius and velocity exist, see Fig. |8j Also, crack branching is predicted for high 
driving forces due to a secondary ATG instability. Since the crack blunts with 
higher driving force and the concentration of stresses in the tip region, unstable 
ATG modes should "fit" into the tip region beyond a critical driving force. At that 
stage, the simulations were not yet quantitative, i.e. the results for this model, 
which is based on sharp interface equations, still depended on the system size and 
the phase field interface thickness. 

Due to the fact that the problem naturally contains several very different length- 
scales, which are the macroscopic system size, the tip scale, the phase field interface 
thickness and the grid spacing for a numerical implementation, the required scale 
separation is difficult to achieve. An alternative sharp interface description for the 
steady state regime was developed in [104] . This approach has the advantage that 
it removes both the phase field interface thickness and the finite system size from 
the description. The idea is that, in contrast to the singular stress field of a sharp 
crack, see Eq. Q, also higher order modes with a stronger divergency at the crack 
tip can be present for a blunt tip: 

a i3 = ( Kf i3 {6) + + Aj-^P- + . . .V (16) 

V2vrr \ r r 1 ) 

Here, the angular dependencies fyk are again known functions, since they are 
eigenmodes of a straight cut, and the weight factors serve as expansion coeffi- 
cients. The higher order modes cannot be excluded, since the usual argument that 
they carry an infinite elastic energy is not valid for a crack with finite crack tip 
radius which serves as a cutoff. Therefore, the stress field becomes a superposition 
of all these modes, and only the slowest decaying one is associated with the stress 
intensity factor, that is related to the driving force. The other expansion coefficients 
Ah are used to solve the linear elastic problem of traction free crack surfaces for 
a steady state crack. Then the shape of the crack and its velocity are determined 
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self-consistently for steady state solutions. The results confirm the basic picture 
of the phase field simulations, i.e. the existence of steady state solutions with a 
velocity appreciably below the sound speed and a finite tip radius. Branching is 
indicated by a transition towards a crack tip with negative curvature in the steady 
state regime. The results differ in detail from [103], and in particular the crack 
velocity is predicted as a partially decaying function of the driving force, which is 
a counterintuitive result. Although the velocity dependence is rather weak, the ap- 
plicability of the model to real materials remains an open issue. Large-scale phase 
field simulations with a careful extrapolation to a full scale separation and infinitely 
large systems confirm these results [105J. 

In the preceding phase field and sharp interface models the crack surface is 
considered as the only source of dissipation, in agreement with the concept of small 
scale yielding. It is however known that for many materials e.g. plastic effects and 
the formation of defects play a crucial role in the vicinity of the crack tip. As 
a first step towards an incorporation of such effects crack growth in viscoelastic 
media has been considered in [106] . Here, the viscous damping coefficients introduce 
together with the surface diffusion coefficient a timescale, which finally provides 
selection for both the crack tip velocity and the tip radius. In contrast to the inertia 
limited regime, the growth velocity is a monotonically increasing function of the 
driving force up to the point of crack branching. Especially for mode I cracks, the 
contribution of the viscous damping turns out to be dominant in comparison to the 
dissipation at the crack front, and leads to a renormalization of the Griffith point; 
below this apparent growth threshold (but above the literal Griffith point from 
energy balance), the crack grows only very slowly. The higher the admixture of a 
mode III loading, the faster the crack propagates in the low driving force regime, 
which may indicate a front instability. First attempts towards linking the inertia 
and bulk damping limited regimes is done in [107J , employing both sharp interface 
and phase field methods. Similarly, a crack tip scale selection can be induced by 
plastic yielding [108J. 

7. Outlook 

Different phase-field approaches for fracture have been developed and analyzed to 
various degrees over the last decade. Crack propagation at low speed is governed es- 
sentially by macroscopic linear elasticity with traction free boundary conditions on 
the fracture surfaces. Therefore simple phenomenological phase field descriptions 
of fracture, such as the KKL model [40] . are generally able to describe this evolu- 
tion quantitatively because the dynamics does not depend sensitively on details of 
the process zone physics. The analysis of such models has reproduced known crack 
propagation laws, such as the principle of local symmetry for isotropic media, and 
has provided a fertile ground for extending those laws to anisotropic media within 
the traditional energetic framework of continuum fracture mechanics. Furthermore, 
numerical examples to date have demonstrated the ability of those models to de- 
scribe quantitatively crack kinking and oscillatory instabilities with biaxial loading 
and during thermal fracture. 

For fast moving inertial cracks, phase field models have been able to reproduce 
crack branching instabilities seen experimentally. However, unlike for slow moving 
cracks, the threshold velocity for branching depends generally on details of the pro- 
cess zone description, which remains largely phenomenological. Therefore, while it 
is encouraging that this threshold is in the observed range of some fraction of the 
wave speed, it is not yet clear how to relate quantitatively phase field model pre- 
dictions of dynamical branching instabilities to experimental observations. Phase 
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field models that incorporate a more realistic description of the process zone are 
presently needed to make predictions more broadly applicable. Inertial cracks have 
also been studied using a different class of phase field models that treat fracture 
as the late stage of the morphological instability of a stressed solid |1U3| . Those 
models have the advantage of having a well-defined sharp- interface limit, which 
makes it possible to derive explicit scaling predictions for the crack velocity. Those 
predictions have highlighted interesting analogies and differences between the tip 
selection problems for cracks and dendritic crystals. However the physical descrip- 
tion of fracture in those models is also too idealized to make quantitative predictions 
for real materials. 

Effects that are related to the discrete nature of the material, including lat- 
tice trapping effects, which produce a velocity gap and also affect macroscopically 
observable properties |109| . are not contained in the present phase field models. 
Incorporating such effects in a continuum formulation presents a major challenge. 
Furthermore, in many materials crack propagation is accompanied by nonlinear 
elastic effects {110H112] . severe plastic deformation, as well as dislocation emission 
that is believed to play a key role in the brittle-to-ductile transition. Phase field ap- 
proaches that incorporate dislocations explicitly and the phase field crystal model 
have the potential to model those processes, but investigations of those models are 
still in very early stages for fracture. Another promising approach to tackle this 
transition is to modify the KKL model to incorporate a continuum description of 
plasticity. 

Finally, phase field models for crack growth have been studied so far only in two 
dimension, with the exception of Ref. |113j . where three-dimensional crack front 
instabilities under mixed mode loading are investigated. The computational cost for 
3D simulations is very high, and eventually grid adaptivity and highly parallelized 
schemes have to be used to provide efficient algorithms. 
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